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ABSTRACT 

It  is  widely  recognized  that  spectral  estimation  serves  as  a  powerful 
digital  signal  processing  tool  in  such  diverse  areas  as  radar,  array  processing, 
adaptive  filtering,  seismology,  and,  speech  processing.  In  this  paper,  a 
novel  method  of  rational  spectral  estimation  shall  be  presented  which  possesses 
a  number  of  admirable  properties:  (1)  it  has  an  elegant  algebraic  structure, 

(ii)  its  spectral  estimation  performance  has  been  empirically  found  to  be 
superior  to  such  contemporary  methods  as  the  maximum  entropy  and  Box-Jenkins 
methods, and  (ill)  it  is  implementable  by  a  "fast  algorithm"  which  is  com¬ 
putationally  competitive  with  recently  developed  fast  LMS  algorithms.  Taken 
in  combination,  these  properties  mark  this  new  method  as  being  a  primary 
spectral  estimation  tool  in  those  challenging  applications  requiring  high  per¬ 
formance  spectral  estimation  in  a  real  time  setting. 

The  new  method  is  predicated  on  selecting  the  parameters  of  a  rational 
spectral  model  so  that  the  underlying  Yule-Walker  equations  are  best  approxi¬ 
mated  over  an  extended  time  span.  This  approach  was  incorporated  by  Cadzow 
in  evolving  an  effective  pole-zero  spectral  model  (e.g.,  see  refs.  [1],[2], 

[22]).  More  recently,  Ogura  and  Yoshida  have  adopted  this  same  procedure  for 
obtaining  a  more  specialized  all  pole  spectral  model  [23]. 


This  work  was  supported  in  part  by  the  Office  of  Naval  R^yareh, 

The  Statistics  and  Probability  Program  under  Contract  N0(Wl4-8O-C-0303/ 


I.  INTRODUCTION 


In  many  signal  processing  applications,  it  is  necessary  to  estimate 
the  power  spectral  density  which  characterizes  a  wide-sense  stationary 
time  series  (x(n)-}.  This  estimate  is  to  be  based  upon  a  finite  set  of 
time  series  observations  which  is  here  taken  to  be  the  n  contiguous 
measurements 

X1 » x2 »  •*•»  XQ  (1) 

A  standard  procedure  for  obtaining  a  spectral  estimate  is  to  first  postu¬ 
late  a  parametric  spectral  model,  and,  then  to  use  the  above  time 
series  observations  for  fixing  the  model's  parameters.  In  light  of  this 
given  observation  set,  it  is  clear  that  the  spectral  model  used  must  be 
dependent  on  only  a  finite  set  of  parameters.  Without  doubt,  the  most 
widely  employed  model  is  the  so-called  rational  spectral  density  model  as 
specified  by 


1-0 +  +  •••  +y-J<1“l2 

|l  +  a1e”^u  +  ...  +  a  e~^ph)|  2 
l  P 


(2) 


This  spectral  density  model  is  generally  referred  to  as  an  autoregressive- 
moving  average  (ARMA)  model  of  order  (p,q).  The  importance  ascribed  to 
this  model  arises  from  the  fact  that  any  continuous  spectral  density  can 
be  approximated  arbitrarily  closely  by  a  rational  function  of  form  (2)  by 


Although  this  ARMA  model  is  the  most  general  rational  spectral 


density  model,  the  more  specialized  autoregressive  (AR)  model  for  which 
q  ■  0,  and,  the  moving  average  (MA)  model  for  which  p  ■  0  have  received 
the  preponderance  of  attention  from  researchers  and  users  alike.  For 
Instance,  the  maximum  entropy,  one-step  predictor,  and,  autoregressive 
methods  have  been  developed  for  efficiently  estimating  an  AR  model's  a^ 
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coefficients.  Similarly,  the  periodogram  approach  and  its  variants  have 


been  developed  for  efficiently  estimating  a  MA  model's  coefficients. 

The  interested  reader  will  find  excellent  treatments  of  these  and  other 
rational  spectral  estimation  methods  in  Haykin  [4]  and  Childers  [5]. 

The  primary  reasons  for  concentrating  on  AR  and  MA  spectral  models 
have  been  that:  (i)  they  have  proven  effective  in  achieving  satisfactory 
spectral  estimation  performance,  and,  (ii)  these  particular  modeling 
approaches  result  in  tractable  analyses  and  efficient  implementations. 

It  is  widely  recognized,  however,  that  an  ARMA  spectral  model  will 
generally  lead  to  an  improved  performance  in  which  fewer  model  para¬ 
meters  are  used.  In  recognition  of  this  fact,  a  variety  of  procedures 
have  been  developed  for  generating  ARMA  models.  These  include  the 
whitening  filter  approach  which  is  typically  iterative  in  nature, 
generally  slow  in  convergence,  and,  usually  requires  a  rather  large 
number  of  time  series  observations  to  be  effective  (e.g.,  see  refs.  [6] 
and  [7]).  More  desirable  closed  form  procedures  which  overcome  most  of 
these  deficiencies  have  been  proposed.  These  include  the  so-called 
Box-Jenkins  method  and  its  variants  [8]-[10],  and,  more  recently,  Cadzow 
has  developed  the  "high  performance" method  [1]  and  [2].  Although  this 
latter  method  has  provided  excellent  spectral  estimation  performance  when 
compared  to  its  ARMA  and  AR  model  counterparts,  its  computational 
efficiency  is  relatively  inferior. 

Recently,  attention  has  been  directed  towards  developing  "fast" 
spectral  estimation  algorithms.  These  algorithms  are  typically  based  on 
the  divide  and  conquer  approach.  As  an  example,  it  is  possible  to  use 
this  approach  for  estimating  the  autoregressive  coefficients  of  a  pc^ 
order  AR  model  with  the  number  of  required  additions  and  multiplications 
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being  on  the  order  of  p  log(p)  (e.g. ,  see  refs.  [ 11 ]— [13 ] > .  These  fest 
algorithms  offer  the  exciting  possibility  of  providing  an  adaptive  spec¬ 
tral  estimator  in  which  the  spectral  model's  parameters  are  updated  as 
new  time  series  observations  become  available.  Unf ortunately ,  the  imple¬ 
mentation  of  these  fast  algorithms  tend  to  be  rather  complex  and  a 
relatively  large  value  for  the  model  order  parameter  p  is  required  before 
the  computational  complexity  p  log(p)  is  approached.  It  is  believed  that 
future  developments  will  alleviate  these  shortcomings. 

With  these  thoughts  in  mind,  it  is  apparent  that  there  is  a  pressing 
need  for  a  high  performance  ARMA  spectral  estimation  method  which  will 
possess  the  computational  efficiency  of  the  fast  algorithms.  In  this 
paper,  a  novel  ARMA  modeling  procedure  for  achieving  these  two  objectives 
will  be  presented.  The  primary  objective  of  this  paper  will  be  that  of 
developing  the  basic  description  of  this  method  and  illustrating  its 
performance  behavior  by  means  of  a  number  of  comparative  examples.  The 
details  of  how  one  may  implement  this  algorithm  so  as  to  attain  the 
aforementioned  fast  computational  complexity  will  be  given  in  a  subse¬ 
quent  paper  [14], 

We  shall  herein  consider  primarily  the  task  of  estimating  the  ARMA 
model's  a^  autoregressive  coefficients.  This  estimation  will  be,  to  a 
large  extent,  motivated  by  the  well  known  Yule-Walker  equations  which 
characterize  stationary  ARMA  time  series.  The  Yule-Walker  equation  develop¬ 
ment  and  other  fundamental  concepts  will  now  be  briefly  reviewed. 


II.  PRELIMINARY  CONCEPTS 

In  chis  section,  a  number  of  concepts  which  are  central  to  the 
development  of  this  paper's  ARMA  model  spectral  estimation  method  are  pre¬ 
sented.  Undoubtably  the  key  notion  is  contained  in  the  Yule-Ualker 
equations  which  characterize  rational  stationary  time-series. 


Yule-Walker  Equations 

It  is  readily  shown  that  the  stationary  random  time  series  whose 
spectral  density  is  specified  by  relationship  (2)  can  be  modeled  as  being 
the  response  of  the  causal  ARMA  system 


P  <1 

\  +  l  Vk-i  *  l  Vk-i  (3) 

i-1  i-o 

to  a  zero  mean  white  noise  excitation  {e  }  whose  individual  terms  have 

n 

variance  one.  The  autocorrelation  description  of  this  system  is  readily 
obtained  by  first  multiplying  each  side  of  expression  (3)  by  the  entity 
x£_m  and  then  taking  the  expected  value.  This  results  in  the  well  known 
Yule-Walker  equations  as  specified  by 


r  (m)  +  l  a.r  (m-i)  -  0 
i-1 


for  m  >  q 


(4) 


In  thisexpression,  the  entries  rx(m)  denote  the  time  series'  auto¬ 
correlation  elements  as  given  by 

rx(m)  -  E{xkx^_m>  (5) 

in  which  *  and  E  denote  the  operations  of  complex  conjugation  and 
expected  value,  respectively. 

In  what  is  to  follow,  the  characteristic  Yule-Walker  equations  (4) 
will  serve  as  the  basis  for  estimating  the  ARMA  model's  a^  and  b^  co¬ 
efficients  from  the  given  set  of  time  series  observations  (1).  The  method 
to  be  described  has  the  appealing  property  of  being  both  computationally 
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efficient  to  implement,  end,  adaptive  in  nature.  Namely,  as  the  new 
time  series  element  x^^  becomes  available,  it  will  be  possible  to 
efficiently  update  the  AFMA  model's  optimal  autoregressive  coefficients 
corresponding  to  the  n  data  length  set  so  as  to  obtain  the  optimal 
autoregressive  coefficients  for  the  new  n+1  data  length  set.  This  capa¬ 
bility  is  essential  if  real  time  spectral  estimation  is  to  be  achieved. 


Down  and  Up- Shift  Operators 

In  the  spectral  estimation  method  to  be  presented,  extensive 
utilization  of  the  down  shift  operator  S  will  be  made.  This  operator 
is  formally  defined  by 

t 

Sx  *  [0,x1,x2,  xn_^^  (6) 


in  which  x  is  an  n*l  vector  as  given  by 

x  -  [xl’x2’  xn3’  (7) 

where  the  prime  symbol  denotes  vector  transposition.  As  its  name  implies, 
operator  S  downshifts  by  one  unit  the  elements  of  the  vector  upon  which 
it  operates  and  inserts  a  zero  into  the  vacated  first  component  position. 
The  downshift  operator  has  the  following  n*n  matrix  representation 


0 

1  0 


o 


(8) 


SO 

L  1  °. 

When  this  operator  is  applied  sequentially  m  times  to  the  vector  x. 


it  is  clear  that  a  down-shift  of  m  units  will  arise  as  follows 


(9) 


m  1 

S  x  *  [0,  0,  ....  0,  x, ,  x~  ....  x  ] 

—  1  * 12  n-m 

m  zeroes 

where  the  integer  m  is  taken  to  be  nonnegative. 

In  a  similar  fashion,  the  up-shift  operator  as  formally  defined  by 
S  ^x  *  [x2,  x3,  ....  xn,  0]  (10) 

will  also  play  a  key  role.  We  have  here  used  an  imprecise  notation  for 
the  up-shift  operator  in  the  sense  that  S  ^ does  not  designate  the  inverse 
of  operator  S.  As  a  matter  of  fact,  the  inverse  of  S  does  not  exist. 

The  primary  reason  for  using  the  notation  S  ^  will  be  due  to  the  notational 
simplification  which  thereby  arises.  Namely,  the  operator  Sm  will  serve 
the  dual  role  of  corresponding  to  a  m  down-shift  operator  when  the 
integer  npO  while  corresponding  to  an  up-shift  operator  whenever  m<0.  In 
using  the  operators  S  and  S  \  it  is  important  to  appreciate  the  fact 
that 

SmSn  »  iff  m  &  n>o  or  m&n<o  (11) 

but  that  when  the  integers  m  and  n  are  of  opposite  sign,  this  equality 
will  be  invalid.  As  a  final  comment  it  is  a  simple  matter  to  show  that 
the  n*n  matrix  S  ^  is  specified  by 

S-1  *  S  (12) 

that  is,  the  up-shift  operator  is  simply  the  matrix  transpose  of  the  down¬ 
shift  operator. 

Displacement  Rank 

The  ability  to  achieve  a  fast  algorithmic  spectral  estimator  will  be 
dependent  upon  one  taking  advantage  of  the  near  Toeplitz  structure  of 
certain  matrices  which  arise  in  the  procedures  to  be  described.  In 
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particular,  the  degree  to  which  the  p*p  matrix  A  is  Toeplitz  in  structure 
is  measured  by  its  "displacement  rank"  [15].  The  displacement  rank  of 
matrix  A  is  formally  given  by 

a(A)  *  min[a_(A),  a+(A)]  (13) 

where 

a_(A)  *  rankfA  -  S  A  S  ]  (14) 

a+(A)  -  rankfA  -  s'a  S]  (15) 

It  is  a  simple  matter  to  show  that  the  displacement  rank  of  a  Toeplitz 
matrix  is  less  than  or  equal  to  two.  Thus  a  matrix  whose  displacement 
rank  is  near  two  is  said  to  be  close  to  Toeplitz  in  structure. 

In  solving  the  system  of  p  linear  equations  in  p  unknowns  as 
specified  by 

Ax  =  y 

one  may  use  standard  matrix  inversion  methods  to  arrive  at  the  solution. 

3  3 

These  routines  will  take  on  the  order  of  p  (i.e.,  <j>(p  ))  multiplication 

and  addition  operations  to  obtain  that  solution.  If  the  displacement 

rank  of  matrix  A  is  a,  however,  one  may  use  a  generalized  version  of  the 

2 

Levinson  algorithm  to  obtain  the  solution  using  <J>(ap  )  computations  [16]. 
If  a  is  sufficiently  smaller  than  p,  it  then  follows  that  a  significant 
computational  savings  can  be  realized  by  adopting  the  generalized  Levinson 
algorithm  solution  procedure.  In  what  is  to  follow,  use  will  be  made  of 


this  observation. 


III.  ARMA  SPECTRAL  ESTIMATION 


In  this  section,  a  novel  procedure  for  estimating  an  ARMA  model's 
autoregressive  coefficients  shall  be  presented.  This  development  is 
begun  by  first  evaluating  the  model  equation  (3)  over  the  integer  set 
p+1  <  k  <  n  to  obtain  the  time  series  relationships 


X  .  , 
p+1 

X 

P 

X  i  • 
p-1 

.  . 

V 

Gp+1 

£  •  •  •  £  . 
p  p-q+i 

b0 

X  .  _ 
p+2 

+ 

p+1 

X  .  . 

p 

.  .  x 2 

a2 

SB 

ep+2 

£p+l  *  ’  ep-q+2 

bl 

x 

n 

X  -i 
n-1 

X  „  . 

n-2 

.  •  X 

n-p 

a 

P 

£ 

n 

n-1  n-q 

b 

q 

(16a) 


In  what  is  to  follow,  it  will  be  convenient  to  represent  these  relation¬ 
ships  in  the  compact  vector  space  format 

x  +  X  a  *  6  _b  (16b) 

where  x,  a  and  _b  are  n-p,  p,  and  q+1  column  vectors,  respectively,  while 
X  and  6  are  (n-p)xp  and  (n-p)x(q+l)"Toeplitz  type" matrices,  respectively. 
The  entries  of  these  individual  vectors  and  matrices  are  directly  obtained 
upon  examination  of  relationships  (16a)  and  (16b). 

It  is  now  desired  to  utilize  relationship  (16b)  in  conjunction  with 
the  Yule-Walker  equations  (4)  to  effect  a  procedure  for  estimating  the  ARMA 
model's  autoregressive  coefficients.  As  we  will  see,  this  objective  is 
attained  by  first  introducing  the  following  (n-p)xt  Toeplitz  type  matrix 
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Y  = 


(17) 


X  X  •  •  •  X  .  ,  <1 

p-q  p-q-1  p-q-t+1 

p-q+1  p-q  p-q-t+2 

•  *  * 

•  •  • 

X  -  X  «  ♦  •  •  X 

n-q-1  n-q-2  n-q-t 

in  which  the  convention  is  adopted  of  setting  to  zero  any  matrix  entry  x^ 
for  which  k  lies  outside  the  observation  index  interval.^  The  integer  t, 
which  specifies  the  number  of  columns  of  matrix  Y,  will  be  found  to 
correspond  to  the  number  of  Yule-Walker  equations  that  are  being  approxi¬ 
mated  (i.e.,  relationship  (4)  for  q<m<_q  +  t).  A  discussion  on  how  one 
selects  the  value  of  t  will  be  shortly  given. 

To  achieve  the  desired  Yule-Walker  equation  approximation,  let  us  now 
premultiply  each  side  of  relationship  (16b)  by  the  complex  conjugate  trans¬ 
pose  of  matrix  Y  as  denoted  by  YA  to  yield 

Yfx  +  Y+Xa  =  Yf€b  (18) 

Upon  taking  the  expected  value  of  this  system  of  t  equations,  it  is  found 
that  for  the  case  when  the  ARMA  model's  numerator  order  is  greater  than 
or  equal  to  its  denominator  order  (i.e.,  q>_p),  this  yields 


A  more  generalized  version  of  this  modeling  process  may  be  obtained  upon 
substituting  the  integer  s  for  p  wherever  it  appears  in  relationship  (16) 
and  (17).  For  simplicity  of  presentation,  the  special  case  s  ■  p  is 
herein  considered. 


In  either  order  case,  it  is  therefore  seen  that  the  system  of  equations 

(18)  provides  an  "unbiased  estimate  of  the  first  t  Yule-Walker  equations. 

It  is  important  to  note  that  the  right  side  vector  Y^Sb  has  an  expected 

value  of  (the  t*l  zero  vector).  This  is  a  direct  consequence  of  the 

ARMA  model's  causality  and  the  whiteness  of  the  excitation  process  which 

renders  E{x  e  *}  5  0  for  all  m  <  k. 
m  k 

With  these  thoughts  in  mind,  an  appealing  procedure  for  selecting 
the  ARMA  model's  p  autoregressive  coefficients  is  suggested.  Namely, 
they  will  be  selected  in  a  manner  so  as  to  best  "satisfy"  the  system  of 
linear  equations 

Y+x  +  YfXa  -  6  (20) 

which  constitutes  the  aforementioned  unbiased  estimate  of  the  first  t 
Yule-Walker  equations.  In  a  real  sense,  this  selection  process  constitutes 
a  procedure  which  is  "most"  consistent  with  the  given  time  series  obser¬ 
vations  and  the  hypothesized  ARMA  model  of  order  (p,q).  The  elements  of  the 
vector  Y*x  ari  matrix  Y+X  for  this  "unmodified"  method  are  given  in  Table  1. 

Relationship  (20)  represents  a  system  of  t  linear  equations  in  the 
p  autoregressive  coefficient  unknowns.  As  such,  a  "best"  solution  to 
this  system  of  equations  is  meaningful  only  in  those  cases  whereby  t  _>  p. 

In  what  is  to  follow,  it  will  be  then  assumed  that  t  _>  p,  and,  that  the 
rank  of  the  t*p  matrix  Y^X  is  p.  Under  these  constraints,  it  is  clear 
chat  the  system  of  equations  (20)  will  be  consistent  when  t  *  p  and  will 
be  generally  inconsistent  when  p  <  t  s  n-q-1.  The  upper  limit  on  t 
(i.a.,.  n-q-i)  constitutes  the  largest  Yule-Walker  equation  approximation 
which  is  consistent  with  the  given  ARMA  model  and  che  data  of  length.  With  ' 
in  mind,  let  us  now  separately  treat  the  cases  t  -  p  and  p  <  t  <_ n-q-1. 
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In  this  case,  Che  p*p  system  of  linear  equations  (20)  has  a  unique 
solution.  If  standard  matrix  inversion  techniques  were  used  to  obtain 


this  solution,  the  required  number  of  multiplication  and  addition  compu- 

3  3 

tations  would  be  on  the  order  of  p  (i.e.  i(p  ).  It  is  possible,  however, 
to  take  advantage  of  the  near  Toeplitz  structure  of  the  pxp  matrix  Y+X 
to  achieve  a  more  efficient  solution  procedure.  Specifically,  it  is 
readily  shown  that  the  displacement  rank  of  this  matrix  is  less  than  or 
equal  to  4.  This  being  the  case,  one  may  then  use  a  generalized  Levinson 
algorithm  to  obtain  the  solution  to  relationship  (20)  using  on  the  order 
of  4p  computations.  Thus,  the  autoregressive  coefficient  selection 
procedure  as  given  by  equation  (20)  is  computationally  competitive  with 
the  maximum  entropy  method.  The  spectral  estimation  performance  of  this 
ARMA  model  procedure,  however,  typically  provides  substantial  improve¬ 
ment  over  the  maximum  entropy  method. 

Relationship  (20)  bears  some  resemblance  to  the  first  iterate  of  the 
Box-Jenkins  method  [8]  of  ARMA  spectral  estimation  in  the  sense  that  the 
first  p  Yule-Walker  equations  (i.e.  equation  (4)  with  q<m<_q+p)  are 
being  approximated.  Upon  closer  examination,  however,  it  is  found  that 
the  standard  Box-Jenkins  method  does  not  possess  the  algebraic  structure 
which  enables  one  to  utilize  sophisticated  LMS  concepts  to  effect  a  super 
efficient  solution  procedure  for  solving  (20).  As  a  final  note,  it  is  a 
simple  matter  to  show  for  the  special  case  q  -  0  and  t  *  p,  that  matrix 
Y  *  X  and  relationship  (20)  reduces  to  the  well-known  covariance  method 
for  obtaining  an  AR  spectral  model  (e.g.,  see  refs.  [17]  and  [13]). 


When  c  >  p,  the  system  of  t  linear  equations  (20)  in  the  p  auto¬ 
regressive  coefficients  is  overdetermined.  As  such,  it  will  not  be 
generally  possible  to  find  a  solution  to  expression  (20).  It  is  possible, 
however,  to  select  an  autoregressive  coefficient  vector  a  which  causes 
this  system  of  equations  to  be  best  approximated  in  the  least  squared  error 
sense.  This  entails  the  Introduction  of  the  following  quadratic  functional 

f (a)  -  [Y*x  +  Y*X  a]* A  [Y+x  +  Y+X  a] 

which  measures  the  degree  to  which  the  system  of  linear  equations  (20)  is 

being  approximated.  The  t*t  diagonal  matrix  A  is  taken  to  be  positive- 
2 

semidefinite.  It  is  a  simple  matter  to  show  that  the  autoregressive  co¬ 
efficient  vector  which  minimizes  this  quadratic  function  satisfies  the 
normal  equations 

XfY  A  TX  a  -  -  XfY  A  Y+x  (21) 

A  solution  to  this  linear  system  of  p  equations  then  provides  the  auto¬ 
regressive  coefficient  estimates  for  the  postulated  ARMA  model.  It  is  of 
interest  to  note  that  this  system  of  linear  equations  reduces  to  the 
Cadzow  High  Performance  method  of  ARMA  spectral  estimation  for  the  particu¬ 
lar  selection  q  -  p  (e.g.,  see  refs.  [1]  &  [2]. 

The  improved  spectral  estimation  performance  achieved  in  using  the 
estimation  method  (21)  over  such  competitors  as  the  3ox-Jenkins  method  is  a 

7 

"The  diagonal  elements  X^  of  the  matrix  A  are  typically  selected  to  be 
nonincreasing  (i.e.,  X^k  2.  -'k+lk+l)  90  as  Co  reflect  an  anticipated  in¬ 
crease  in  the  variance  of  the  Yule-Walker  equation  approximations  for 
increasing  k  as  evident  from  relationships  (19).  Choices  of  * (n-q-k) 
‘kk  *  COk  have  been  found  to  yield  satisfactory  performance. 

I  . 


direct  consequence  of  selecting  integer  t  to  be  larger  than  the  minimal 
value  p.  With  the  corresponding  larger  set  of  Yule-Walker  equations  that 
are  being  approximated,  it  intuitively  follows  that  the  modeling  auto¬ 
regressive  coefficients  will  be  less  sensitive  to  Inaccuracies  of  the 
correlation  estimates  as  embodied  in  Y^X  and  Y^x  than  in  the  minimal  value 
situation  when  t  *  p.  This  anticipation  of  a  spectral  estimation  per¬ 
formance  improvement  has  in  fact  been  demonstrated  on  a  rather  large 
number  of  numerical  examples  carried  out  by  the  authors  and  others. 

From  a  computational  viewpoint,  the  primary  advantage  accrued  in 
using  the  herein  developed  autoregressive  coefficient  estimation  procedure 
resides  in  the  specific  algebraic  structure  thereby  obtained.  This 
structure  may  be  used  to  effect  particular  "fast"  algorithmic  solution 
procedures  for  the  autoregressive  coefficient  estimates.  These  fast 
algorithms  are  computationally  much  more  efficient  than  the  aforementioned 
generalized  Levinson  algorithm.  In  order  to  achieve  the  required  structure 
for  these  fast  algorithms,  it  will  be  necessary  to  make  minor  modifications 
on  the  constituent  vector  x  and  matrices  X  and  Y.  In  the  next  section, 
a  discussion  of  some  candidate  modifications  will  be  made. 
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IV.  PRE  AND  POSTMODIFICATION  METHODS 


It  is  possible  Co  realize  significant  computational  savings  in 
the  proposed  ARMA  spectral  estimation  method  over  that  provided  by  the 
generalized  Levinson  algorithm.  This  improvement  will  entail  a  slight 
modification  in  Che  constituent  vector  x  and  matrices  X  and  Y  which 
constitute  the  Yule-Walker  equation  estimations  (20).  Although  the 
suggested  modifications  will  typically  result  in  biased  estimates  of  the 
Yule-Walker  equations,  it  is  shown  that  when  the  data  length  n 
adequately  exceeds  the  order  parameters  p  and  q  (i.e.,  n  >>  p&q) 
that  these  estimates  are  virtually  unbiased. 

In  addition  to  these  vector  and  matrix  modifications,  it  will  be 
necessary  to  restrict  the  integer  t  to  its  minimal  value  of  p  in 
order  to  achieve  the  fast  algorithm  solution  capability.  Unfortunately, 
the  restriction  t  *  p  will  generally  result  in  an  associated  decrease 
in  spectral  estimation  performance.  Thus,  in  obtaining  a  computationally 
fast  algorithmic  solution  procedure  for  the  a^  coefficients,  an  accom¬ 
panying  sacrifice  in  spectral  estimation  performance  is  the  price 
being  paid.  One  must  therefore  carefully  consider  the  ramifications  of 
this  tradeoff  for  any  given  application.  Fortunately,  the  degradation  in 
performance  is  not  great  for  many  relevant  applications. 

With  these  remarks  in  mind,  we  shall  now  consider  the  aforementioned 

modifications  in  the  constituent  x,  X,  and  Y  entries.  Although  a 

~  / 

t 

countless  number  of  modifications  are  possible,  we  will  be  concerned 
exclusively  with  the  so-called  "premodification"  and  "postmodification" 
procedures.  Additional  modification  procedures  are  now  being  examined  and 
will  subsequently  be  reported  upon. 
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(i)  Pr«modificatioo  Method 

I 

In  Che  premodification  method,  the  vector  x,  and,  matrices  X  and  Y 

used  in  expressions  (16)  and  (17)  are  modified  as  follows.  The  vector  x 

becomes  the  given  observed  time-series  vector 

» 

■  [Xj,  x2,  ...  xnl  (22a) 

while  the  n*p  matrix  X  has  the  p  column  vector  structure 

Xx  -  [Sxx  i  s2^  :  •  •  •  :  s^)  (22b) 

and  the  n*t  matrix  Y  has  the  t  column  vector  structure 

Y!  -  [sq+1Xi  :  sq+2x1  :  •  •  •  :  sq+tx1]  (22c) 

where  S  is  the  aforementioned  n*n  down-shift  matrix.  When  these 
entries  are  substituted  into  relationship  (20),  the  so-called  premodifi¬ 
cation  method  of  estimating  the  ARMA  model’s  autoregressive  coefficients 
is  at  hand.  This  method  is  suggestively  referred  to  as  being  premodified 
since,  among  other  things,  the  original  system  of  equations  (16)  has 
appended  to  it  a  set  of  p  new  equations  as  represented  by  the  first  p 
rows  of  matrix  X^  (the  last  n-p  rows  of  matrix  X^  are  equal  to  matrix  X). 

These  first  p  rows  form  a  lower  triangular  matrix  and,  as  such,  an  un-. 
desirable  "transient"  effect  on  the  ARMA  modeling  is  introduced.  It  is  pre¬ 
cisely  because  of  this  triangular  structure,  however,  that  we  are  able  to 
obtain  a  fast  computational  algorithm  for  achieving  estimates  of  the  ARMA 
models'  autoregressive  coefficients[14] . 

If  the  modifications  as  embodied  in  expressions  (22)  are  substituted 
into  relationship  (20),  the  following  new  approximation  of  the  first  t 
Yule-Walker  equations  is  obtained 

Y1*X1  -  +  V-l  “  9  (23) 
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Upon  caking  Che  expecced  value  of  chis  relacionship,  ic  is  found  ChaC  for 
che  cypical  ARMA  model  order  selecCion  case  q  <_  p,  che  following  biased 
escimaces  of  che  Yule-Walker  equacions  resulc 
m  P 

(n-m)  l  a.r  (m-k)  +  £  (n-k)a,r  (m-k)  -0  q<m<p  (24a) 

k-0  K  x  k-nrfl  x 

P 

(n-m)  £  a.r  (m-k)  ■  0  P£®Ji<l+t  (24b) 

k-0  K  x 

in  which  ag  -  1.  Ic  is  Co  be  noced  ChaC  Chese  escimaces  are  virCually 
unbiased  whenever  n  >>  p.  In  a  similar  manner,  for  che  case  q  >  p,  che 
expeccacion  operacion  resulCs  in  che  unbiased  resulCs 
P 

(n-m)  l  a  r  (m-k)  -  0  q<m£q+t  (24c) 

k-0  K  x 

Thus,  in  any  case,  Che  linear  equacion  (23)  provides  a  suiCably  good 
approximacion  Co  Che  underlying  Yule-Walker  equacions. 

Special  Case  :  C  -  p 


The  primary  advancage  accrued  by  inCroducing  che  premodified  modifi- 

cacion  is  compuCaCional  in  nacure.  IC  is  a  simple  maCCer  Co  show  ChaC, 

in  Che  special  case  C  -  p,  che  pxp  macrix  Y^X^  has  a  displacemenC 

rank  less  Chan  or  equal  Co  3.  As  such,  a  generalized  Levinson  algorichm 

2 

may  be  developed  for  solving  relacionship  (23)  which  requires  <p  (3p  )  compuCacioni 
More  imporcancly,  however,  is  che  face  ChaC  che  parcicular  scruccure  of 


maCrices  and  Y^  will  enable  one  Co  uCilize  a  fasC  algorithmic  pro¬ 

cedure  Co  obcain  chis  solucion  using  $  (p)  compuCations.  This  algorichm 
is  described  in  reference  [14]. 

(ii)  Postmodification  Mechod 

When  employing  che  posCmodif icaCion  mechod,  che  column  vecCor  x  is 


replaced  by  the  n*l  vector 
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*2  "  S_P— ]_  (25a) 

where  x^  is  given  by  relationship  (22a).  The  matrix  X  is  replaced  by 
the  n*p  matrix 

^2  *  [S  ^x^  S  p+2x^I  •  •  *  ‘  x^]  (25b) 

while  the  matrix  Y  is  replaced  by  the  nxt  matrix 

y2  -  [s‘p+q+1x1:  s_p+q+2x1:  •  •  .  :  s'^+\)  (25c) 

Upon  substitution  of  these  entries  into  expression  (20),  the  so-called 
postmodification  method  of  estimating  the  ARMA  model's  autoregressive  co¬ 
efficients  is  obtained.  This  method  is  referred  to  as  being  postmodified 
since,  among  other  things,  the  original  system  of  equations  (16)  are  here 
being  supplemented  with  p  new  equations  as  represented  by  the  last  p 
rows  of  matrix  X2  (the  first  n-p  rows  of  matrix  X2  are  equal  to  matrix 
X) .  These  last  p  rows  form  a  lower  triangular  matrix  and  introduce  an 
undesirable  transient  effect  on  the  resultant  autoregressive  coefficient 
solution.  The  primary  advantage  obtained  in  using  the  postmodification 
method  is  computational  in  nature  as  will  now  be  briefly  described. 

With  these  postmodified  matrix  and  vector  entries,  the  following 
approximation  to  the  first  t  Yule-Walker  equations  is  obtained 

Y2+X2  a  +  Y2+x2  -  9  (26) 

The  degree  to  which  this  linear  system  of  equations  approximates  the  Yule- 
Walker  equation  may  be  measured  by  taking  its  expected  value.  This  is 
found  to  yield  for  the  more  typical  order  selection  case  q  <_  p 
m  P 

l  (n-p+k)a.r  (m-k)  +  (n-p-hn)  £  a,r  (m-k)  -  0  q<m<^p  (27a) 
k-0  X  k-mfl  K  x 

P 

l  (n-nrt-k)a.r  (ra-k)  -  0 
k-0 
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p<m<_q+t  (27b) 


while  for  che  less  typical  case  q>p,  one  finds  that 
P 

l  (n-nrt-k)a.  r  (m-k)  -  0  q<m£q+t  (27c) 

k*0  K  X 

Thus,  relationship  (26)  represents  a  biased  estimate  of  the  Yule-Walker 
equations.  It  is  to  be  noted,  however,  that  when  the  number  of  time 
series  observations  n  adequately  exceeds  the  ARMA  model  parameter  p, 
relationship  (26)  provides  an  almost  unbiased  estimate. 


Special  Case  :  t  *  p 


when  t  -  p,  it  is  a  simple  matter  to  show  that  the  displacement 
r.T>rl-  of  the  pxp  matrix  Y2+X2  is  ^ess  than  or  equal  to  three.  A 
ger  r’  »\ized  Levinson  algorithm  may  therefore  be  used  to  solve  relation- 

2 

snip  (26)  which  requires  <p ( 3p  )  computations.  It  is  possible,  however, 
tv  utilize  the  algebraic  structure  of  relationship  (26)  to  generate  a 
fast  algorithmic  solution  procedure  that  requires  $>(plogp)  computations 
as  presented  in  reference  (14]. 

(iii)  Pre-and-Postmodif ication 

If  the  pre  and  postmodification  procedures  are  simultaneously  in¬ 
corporated,  a  third  procedure  arises.  Specifically,  the  column  vector  x 
is  replaced  by  the  (n+p)*l  vector 


[x.  ,  x«,  *  *  ,  x  ,  0, 


,  0] 


(28a) 


p  zeroes 

The  matrix  X  is  next  replaced  by  the  (n+p)xp  matrix 

x3  “  (sx3  :  s2x3  :  •  •  •  :  spx3) 


(28b) 


while  the  matrix  Y  is  replaced  by  the  (n+p)*t  matrix 
Y  *  [s^+^X  *  •  .  •  .  '  g*l+tx  J 


(28c) 
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If  these  entries  are  substituted  into  relationship  (20),  the  so-called 
"Pre-and-Postmodif ication"  method  of  ARMA  model  autoregressive  co¬ 
efficient  estimation  procedure  is  generated.  In  using  this  particular 
modification,  it  is  readily  shown  that  the  original  system  of  equations 
(16)  are  being  supplemented  with  2p  new  equations  as  represented  by 
the  first  and  last  p  rows  of  matrix  (the  centermost  n-p  rows 

of  matrix  are  equal  to  matrix  X).  The  initial  and  final  p  rows 

introduce  an  associated  deleterious  transient  effect.  The  resultant 
algebraic  structure  for  the  system  of  linear  equations  (11),  however, 
result  in  a  significant  computational  savings  in  determining  the  auto¬ 
regressive  coefficients. 

In  using  the  Pre-and-Postmodification  method,  the  approximation  to 
the  first  t  Yule-Walker  equations  is  provided  by  the  linear  system  of 
equations 

Y3+  X3  a  +  Y3fx3  *  0  (29) 

This  system  of  equations  provides  a  biased  estimate  of  the  Yule-Walker 
equations  as  is  made  apparent  upon  taking  its  expected  value.  For  the 
typical  ARMA  modeling  order  selection  of  q  _<  p,  this  expectation  is 
found  to  yield 

m  P 

£  (n-m+k)a,r  (m-k)  +  £  (n+m-k)a,r  (m-k)  =  0  q<m<p  (30a) 

k-o  K  x  k-mfl  X  X 

P 

£  (n-nH-k)a.r  (m-k)  =  0  p<m<^q+t  (30b) 

k*c  x 

*o 

while  for  the  selection  q>p,  we  have 
P 

£  (n-nH-k)a^rx(m-k)  *  0  q  <  m  <  q+t  (30c) 
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Although  the  biased  nature  of  estimate  (29)  Is  revealed  by  expressions 
(30),  it  is  clear  that  when  the  number  of  time  series  observations  n 
adequately  exceed  the  AKMA  model  order  parameter  p,  the  estimate 
becomes  almost  unbiased. 

Special  Case  :  t  *  p 

In  this  case,  the  displacement  rank  of  the  p*p  matrix  is 

found  to  be  less  than  or  equal  to  two.  A  generalized  Levinson  algorithm 

may  then  be  utilized  to  obtain  the  solution  to  expression  (29)  using 
2 

$(2p  )  computations.  As  in  the  other  modified  methods,  however,  one  may 
utilize  even  more  efficient  solution  procedures  to  obtain  this  solution 
using  <J>(plogp)  computations  [14]. 
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The  following  Table  is  useful  when  programming  the  autoregressive 
coefficient  selection  procedure  as  embodied  in  expressions  (20  or  (21) . 


Method 

A  -  Y+X 

aij  f°r  1  ±  i  £  t 

1  1  j  1  P 

b  =  Y+x 

b.^  for  1  <_  i  <_  t 

Unmodified 

aij  *  jJ-^P+k-j  Xp-q+k-i 

n-p 

b .  *  J  x  . .  x  , ,  . 

1  k-1  p+k  P-<l+k~i 

Premodifiad 

'  jj’VjVq-t 

Postmodified 

n 

aij  “  k^1Xp+k-j Xp-q+k-i 

b.  *  y  x  ..  . 

1  k=i  p_^+k_i 

Pre  &  Postmodified 

n+p 

aij  “  J1Xk-jXk-q-i 

bi  "  ^Vk-q-i 

TABLE  1:  Entries  of  the  matrix  Y  X  and  vector  Y+x  for  the  various 

methods.  The  convention  is  here  adopted  of  setting  x^  *  0 
whenever  k  is  not  in  the  observation  index  interval  1  <  k  <  n. 
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Method 

Displacement 

Rank 

Fast  Algorithm 

Computational 

Measure 

Unmodified 

4 

Premodified 

3 

4>(p) 

Post  modified 

3 

<p(p  log  p) 

Pre  &  Postmodified 

2 

<Hp  log  p) 

TABLE  2:  Computational  characterization  of  various  methods  for  the 
minimal  selection  t  =  p. 
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V.  MODEL  ORDER  SELECTION 


An  important  consideration  in  any  rational  spectral  estimation 
method  is  that  of  the  determination  of  the  model's  order.  A  particularly 
attractive  procedure  exists  for  this  objective  relative  to  the  method 
developed  in  this  paper.  It  is  predicated  on  the  observation  that  the 
ARMA  model's  autoregressive  coefficients  are  obtained  by  solving  a 
system  of  p  consistent  linear  equations  in  the  p  autoregressive  co¬ 
efficient  unknowns.  The  form  of  these  equations  is 

A  a  3  b  (31) 

in  which  A  is  an  approximation  of  an  autocorrelation  matrix.  This  is 
exemplified  by  equation  (20)  with  t  =  p  where  A  *  Y^X.  If  the  time 
series  being  analyzed  is  ARMA  of  order  p,  the  matrix  A  conceptually 
becomes  ill-conditioned  when  the  selected  model  order  p  exceeds  the 
time  series  order  p.  Thus,  the  model  order  determination  can  be  achieved 
by  examining  the  conditioning  of  the  matrix  A  as  a  function  of  p.  As 
p  is  increased,  the  appropriate  choice  will  be  that  value  p  for  which 
there  is  a  precipitate  decrease  in  conditioning  when  p  =  p+1.  This 
approach,  as  applied  to  the  high  performance  method  of  ARMA  spectral 
estimation  [1],  was  suggested  by  Pao  and  Lee  [20]. 

There  exists  a  variety  of  procedures  for  measuring  the  conditioning 
of  a  p*p  matrix  A.  One  of  the  more  popular  such  measures  is  provided 
by  the  normalized  determinant  as  defined  by 

n  p  z 

C (A)  ■  det(A)/v  £  l  j a . . j  (32) 

i=l  j«l  J 

where  det(A)  denotes  the  determinant  of  matrix  A.  It  is  to  be  noted  that 
this  normalized  determinant  equals  zero  when  the  rank  of  the  pxp  matrix 
is  less  than  p.  This  method  has  been  found  to  provide  an  effective  model 
order  determinate [20] . 
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VI.  NUMERATOR  DYNAMICS 

A  variety  of  procedures  exist  for  determining  the  numerator 
dynamics  of  an  ARMA  time  series  once  the  AR  coefficients  have  been 
estimated.  In  this  section,  a  procedure  which  has  been  found  to  be 
particularly  effective  shall  be  presented.  It  makes  use  of  the 
governing  ARMA  relationship  (3)  that  models  the  underlying  time 
series. 

In  this  approach  to  estimating  the  numerator  dynamics,  we  first 
introduce  the  so-called  causal  image  of  a  time  series'  autocorrelation 
sequence  as  specified  by 

r*(n)  m  ~  1  rx(°)<S(n)  +  rx(n)u(n)  (33) 

in  which  6 (n)  and  u(n)  designate  the  standard  unit-sample  and  unit- 
step  sequences,  respectively.  Making  use  of  the  complex  conjugate  sym¬ 
metrical  property  of  stationary  autocorrelation  sequences,  it  then  follows 
that  the  autocorrelation  sequence  can  be  uniquely  recovered  from  its 
causal  image  according  to  the  simple  relationship 

r*(n)  =  r*(n)  +  r*(-n)  0*0 

Upon  taking  the  discrete-Fourier  transform  of  tnis  relationship,  it 
follows  that  the  time  series'  spectral  density  is  given  by 

-U  I  * 

Sx(uj)  =  Sx(oj)  +  Sx(w; 

-  2  Re  {  Sx(u)}  (35) 

where  S^Cw)  designates  the  discrete-Fourier  transform  of  the  causal 
image  sequence  (r^Xn)}.  According  to  this  relationship,  one  may  obtain 
the  required  spectral  density  estimate  by  alternatively  estimating  Sx(w). 
This  will  be  the  approach  now  taken. 
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It  is  possible  to  obtain  a  convenient  closed  form  expression  for 


the  function  S^(u).  This  first  entails  introduction  of  an  auxiliary 
sequence  {c^}  whose  elements  are  specified  by 

P 

Cm  *  rx(m)  +  1iiakrx(m_k)  (36) 

k»l 

After  examination  of  the  Yule-Walker  equations  (4)  which  govern  the  ARMA 
time  series  and  noting  the  causality  of  the  sequence  {r^(k)}»  it 

follows  that  this  auxiliary  sequence  is  itself  causal  and  finite  in 
length.  Namely, 

cm  =  0  for  0  >  m  >  u  =•  max(q,p)  (37) 

Upon  taking  the  discrete  Fourier  transform  of  relationship  (36)  and  in¬ 
corporating  identity  (37),  we  have  after  rearrangement  of  terms 


c  +  c,e 
o  1 


-rue 

U 


1  +  a  e  -$-  •  *  •  +  a  e 
1  P 


-jpw 


(38) 


where  it  is  recalled  that  u  *  max(q,p).  If  this  expression  is  then  sub¬ 
stituted  into  relationship  (35),  the  required  formulation  of  the  spectral 
density  is  comploi ad. 

In  summary,  the  procedure  here  presented  for  ascertaining  the 
numerator  dynamics  effect  on  the  overall  spectral  estimate  entails  the 
following  four  step  procedure.  One  first  computes  estimates  of  the  auto¬ 
correlation  elements  rx(0),  rx(l),  *  •  •,  rx(u)  from  the  given  time 
series  observations  (1).  One  may  use  any  procedure  such  as  the  standard 
unbiased  or  biased  autocorrelation  estimate  methods  for  this  task.  Next, 
the  causal  image  elements  r+(k)  for  0  _<  k  <_  u  are  generated  according 
to  relationship  (33).  The  c^  coefficient  estimates  are  then  determined 
from  recursion  (36)  in  which  the  a^  coefficients  generated  from  solving 
(or  approximating)  relationship  (20)  are  utilized.  Finally,  the  resultant 


-23- 


*• 


V 


close  co  zero  for  tn  >  max(q,p),  this  is  indicative  of  an  incompatability. 
This  incompatability  can  be  caused  from  such  factors  as  an  inadequate 
model  order  selection  (p,q),  poor  estimates  of  the  ARMA  model's  a^ 
autoregressive  coefficients,  or,  inaccurate  estimates  of  the  autocorrela¬ 
tion  elements  r  (o) ,  r  (1) ,  •  •  •  ,  r  (y) . 

XX  X 
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VII.  NUMERICAL  EXAMPLES 


ARMA  Model  Example 

The  unmodified  ARMA  modeling  method  of  spectral  estimation,  as  pre¬ 
sented  in  Section  III,  has  been  found  to  possess  a  significantly  superior 
performance  when  compared  to  such  contemporary  alternate  methods  as  the 
periodogram,  maximum  entropy,  and,  the  Box- Jenkins  method  when  applied 
to  'narrow"band  time  series  (i.e.,  sum  of  sinusoids  and  white  noise  [1], 
[2],  [20]).  With  this  in  mind,  the  effectiveness  of  both  the  unmodified 
and  modified  ARMA  modeling  procedures  will  now  be  examined  for  a 

It  tl 

moderately  wide  band  time  series.  In  particular,  we  shall  treat  the 
time  series  as  recently  considered  by  Bruzzone  and  Kaveh  [21].  Specifi¬ 
cally,  their  ARMA  time  series  of  order  (4,4)  is  characterized  by 

xk  -  x£  +  x£  +  0.5ek  (39a) 

1  2 

where  the  individual  time  series  x^  and  x^  are  generated 
according  to 

x^  *  0.4x£_^  -  0.93xk_2  +  £k  (39b) 

4  *  -°-5xLl  “  °'  93xk-2  +  4 
1  2 

in  which  the  ,  and  are  uncorrelated  Gaussian  random 

variables  with  zero  mean  and  unit  variance.  It  then  follows  that  the 
spectral  density  characterizing  the  time  series  (39)  is  given  by 
i  ,-2  .  ,-2 

S  (u>)  »  1 1  -  0. 4e~^u  +  0 . 93e~^  2ai  |  +  1 1  +  0.  Se-^  +  0. 93e_:32a)  |  +  0. 25 

(40) 

Using  the  time  series  description  (39) ,  twenty  different  sampled 
sequences  each  of  length  64  were  generated.  These  twenty  observation 
sets  were  then  used  to  test  various  spectral  estimation  methods.  In 
Figure  1,  the  twenty  superimposed  plots  of  the  ARMA  model  spectral 


estimates  of  order  (4,4)  obtained  using  the  first  iterate  of  the  Box-Jenkins 

k— 1 

method,  and,  this  paper's  unmodified  method  with  X^  «  (0. 95)  and  selections  of 
t  ■  4,  8,  and,  20  are  shown.  For  comparison  purposes,  the  ideal  spectrum 
(40)  is  also  shown.  From  these  plots,  two  observations  may  be  made: 

(i)  the  unmodified  method  with  t  *  4  yields  a  marginally  better  spectral 
estimate  than  the  Box-Jenkins  method,  and,  (ii)  the  unmodified  spectral  esti¬ 
mates  improves  significantly  as  t  is  increased  from  the  minimal  value  4.  This 
latter  observation  is  most  noteworthy  and  indicates  that  the  incorporation 
of  more  than  the  minimal  number  of  Yule-Walker  equations  for  determining 
the  ARMA  model's  autoregressive  coefficients  has  the  anticipated  effect 
of  significantly  improving  spectral  estimation  performance. 

To  test  the  effect  of  an  improper  model  order  selection,  the  same 
twenty  observation  sets  were  next  used  to  obtain  an  ARMA  model  of  order 
(6,6).  The  spectral  density  estimates  which  resulted,  when  using  the 
Box-Jenkins,  and  this  paper's  unmodified  method  with  t  »  6,  12,  and,  30 
are  shown  in  Figure  2.  The  quality  of  the  spectral  estimates  has 
improved  over  that  of  the  proper  order  selection  (4,4)  for  all  methods. 

It  is  noteworthy,  however,  that  the  averaged  normalized  determinant  (32) 

-3 

for  the  unmodified  methods  fell  by  a  factor  of  10  as  p  is  increased 

i 

i 

from  4  to  6  thereby  properly  indicating  an  excessive  model  order  selection.  ] 

Inexplicably,  the  mean  condition  number  for  the  Box-Jenkins  method 
increased  when  p  was  increased  from  4  to  6. 

Next,  the  modification  methods  developed  in  Section  IV  were  applied 
to  these  twenty  different  sampled  sequences  of  length  64  to  obtain  ARMA 
model  spectral  estimates  of  order  (4,4).  The  resultant  spectra  are 
shown  in  Figure  3  where  it  is  apparent  that  only  "a  modest"  degradation  in 
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spectral  estimation  performance  has  accrued  due  to  the  transient  effects 
introduced  by  the  modified  methods.  This  is  indeed  welcomed  news  given 
the  ability  to  implement  these  modified  methods  with  exceptionally  fast 
algorithms.  It  is  to  be  noted  that  the  "postmodified",  and  the  "pre  and 
postmodified"  methods  are  identical  in  this  example. 

As  a  final  example,  twenty  different  sampled  sequences  each  of 
length  200  were  generated  according  to  expression  (39).  With  this  longer 
data  length,  it  was  anticipated  that  an  improvement  in  spectral  estima¬ 
tion  performance  would  result.  A  marked  improvement  is  in  fact  realized 
as  is  made  evident  from  Figure  4  where  the  ARMA  model  spectral  estimates 
of  order  (4,4)  are  shown  for  the  Box-Jenkins  method  and  the  unmodified 
method  for  selections  of  t  ■  4,  S,  and  20. 

AR  Model  Example 

In  order  to  further  demonstrate  the  effectiveness  of  selecting  t 
larger  than  the  minimum  value  p,  let  us  consider  the  time  series  composed 
of  two  sinusoids  at  normalized  frequencies  0.4  and  0.5  in  additive  white 
noise  such  that  the  individual  SNR's  are  lOdB  and  OdB,  respectively,  that  is 
x(n)  ■  /l0  sin(0.4un)  +  sin(0.5iln)  +  c(n) 

Twenty  different  samples  of  this  time  series  each  of  length  64  where  then 
generated.  Next  an  AR  model  of  order  eight  (i.e.,  p»0,  q  *  8)  using  equation 
(21)  for  each  of  these  twenty  samples  where  generated  for  selections  of  t*8, 
and,  t*24.  The  resultant  spectral  estimates  are  shown  in  Figure  5  where  it 
is  clear  that  the  standard  covariance  method  (i.e.,  t  *  8)  cannot  resolve 
the  two  sinusoids  while  the  AR  model  for  which  t  *  24  easily  achieves  the 
desired  resoltuion.  This  is  a  most  noteworthy  observation  in  that  most  AR 
modeling  techniques  (e.g.,  the  covariance,  the  maximum  entropy,  and,  many 
fast  LMA  methods)  either  implicitly  or  explicitly  specify  t«q  and  thereby 
have  an  inherent  inferior  spectral  estimation  performance  basis. 
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VIII.  CONCLUSION 


A  fast  algorithmic  method  for  achieving  an  ABMA  model  of  a  stationary 
random  time  series  has  been  presented.  This  method's  development  is  pre¬ 
dicated  on  the  concept  of  approximating  the  underlying  model's  Yule-Walker 
equations  through  specially  structured  data  dependent  matrices  and  vectors. 
The  chosen  structures  enables  one  to  solve  for  the  ABMA  model's  auto¬ 
regressive  coefficients  using  an  exceptionally  efficient  algorithm.  In 
particular,  the  computational  efficiency  of  this  algorithm  is  competitive 
with  such  AR  model  procedures  as  the  maximum  entropy  and  recently  developed 
fast  LMS  algorithms  [11] -[13]. 

In  addition  to  its  fast  algorithmic  Implementation  capability,  this 
new  method  possesses  a  superior  spectral  estimation  performance.  It  was 
shown  to  be  an  extension  of  the  high  performance  method  of  Cadzow  [1]  and 
[2]  which  has  been  empirically  found  to  provide  superior  spectral  estimates 
when  compared  to  such  alternatives  as  the  maximum  entropy  and  Box-Jenkins 
methods.  Given  its  dual  capability  of  high  spectral  estimation  perform¬ 
ance  and  fast  computational  implementation,  this  new  method  promises  to  be 
a  primary  spectral  estimation  tool. 
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Figure  2:  ARMA  Spectral  Estimates  cf  Order  (5,5) 
Dcta  Length  of  54,  and,  •  -  0.95. 
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Figure  4:  ARMA  Soectral  Estimates  of  Order  (4,4), 
Data  Length  of  200,  and,  \  =  0.35. 


NORMALIZED  FREQUENCE 


FIGURE  5:  Two  Sinusoids  in  Additive  White  Noise  with 
Frequencies  0,5  (10  dB)  ond  0,4  (0  dB). 

AR  Models  of  order  3  with  a  Dote  Length  of  64, 


